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ABSTRACT 



We give a self-contained modern linear stability analysis of a system of n 



^ \ equal mass bodies in circular orbit about a single more massive body. Starting 

O ' with the mathematical description of the dynamics of the system, we form the 



linear approximation, compute all of the eigenvalues of the linear stability matrix, 
and finally derive inequalities that guarantee that none of these eigenvalues have 
Ph| positive real part. In the end, we rederive the result that J.C Maxwell found 

O ' for large n in his seminal paper on the nature and stability of Saturn's rings, 

^ \ which was published 150 years ago. In addition, we identify the exact matrix 

' that defines the linearized system even when n is not large. This matrix is then 

investigated numerically (by computer) to find stability inequalities. Further- 
^ ' more, using properties of circulant matrices, the eigenvalues of the large 4n x 4n 

■ matrix can be computed by solving n quartic equations, which further facilitates 

the investigation of stability. Finally, we have implemented an n-body simulator 
and we verify that the threshold mass ratios that we derived mathematically or 
numerically do indeed identify the threshold between stability and instability. 
Throughout the paper we consider only the planar ra-body problem so that the 
analysis can be carried out purely in complex notation, which makes the equa- 
tions and derivations more compact, more elegant and therefore, we hope, more 
transparent. The result is a fresh analysis that shows that these systems are 
always unstable for 2 < n < 6 and for n > 6 they are stable provided that the 
central mass is massive enough. We give an explicit formula for this mass-ratio 
threshold. 



Subject headings: planets: rings 
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1. Introduction 

One hundred and fifty years ago, Maxwell (1859) was awarded the prestigious Adam's 
prize for a seminal paper on the stability of Saturn's rings. At that time, neither the 
structure nor the composition of the rings was known. Hence, Maxwell considered various 
scenarios such as the possibility that the rings were solid or liquid annuli or a myriad of small 
boulders. As a key part of this last possibihty, Maxwell studied the case of n equal-mass 
bodies orbiting Saturn at a common radius and uniformly distributed about a circle of this 
radius. He concluded that, for large n, the ring ought to be stable provided that the following 
inequality is satisfied: 

mass(Rings) < 2.298mass(Saturn)/n^. 

The mathematical analysis that leads to this result has been scrutenized, validated, and 
generalized by a number of mathematicians over the years. 

We summarize briefly some of the key historical developments. Tisserand (1889) derived 
the same stability criterion using an analysis where he assumed that the ring has no effect on 
Saturn and that the highest vibration mode of the system controls stability. More recently, 
Willerding (1986) used the theory of density waves to show that Maxwell's results are correct 
in the limit as n goes to infinity. Pendse (1935) reformulated the stability problem so that it 
takes into account the effect of the rings on the central body. He proved that, for n < 6, the 
system is unconditionally unstable. Inspired by this work, Salo and Yoder (1988) studied 
coorbital formations of n sateUites for small values of n where the sateUites are not distributed 
uniformly around the central body. They showed that there are some stable asymmetric 
formations (such as the well-known case of a pair of ring bodies in L4/L5 position relative 
to each other — i.e., one leading the other by 60deg). Finally, Scheeres and Vinh (1991) 
extended the analysis of Pendse to find the stability criterion as a function of the number of 
satellites when n is small. The resulting threshold depends on n but for n > 7, it deviates 
only a small amount from the asymptotically derived value. 

In this paper, we give a self-contained modern linear stability analysis of a system of 
equal mass bodies in circular orbit about a single more massive body. We start with the 
mathematical description of the dynamics of the system. We then form the linear approxi- 
mation, compute all of the eigenvalues of the matrix defining the linear approximation, and 
finally we derive inequalities that guarantee that none of these eigenvalues have positive real 
part. In the end, we get exactly the same result that Maxwell found for large n. But, in addi- 
tion, we identify the exact matrix that defines the linearized system even when n is not large. 
This matrix can then be investigated numerically to find stability inequalities even in cases 
where n is not large. Furthermore, using properties of circulant matrices, the eigenvalues 
of the large An x An matrix can be computed by solving n quartic equations, which further 
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facilitates the investigation. Finally, we have implemented an n-body simulator based on a 
leap-frog integrator (see Saha and Tremaine (1994); Hut et al. (1995)) and we verify that 
the threshold mass ratios that we derived mathematically or numerically do indeed identify 
the threshold between stability and instability. 

Throughout the paper we consider only the planar n-body problem. That is, we ignore 
any instabilities that might arise due to out-of-planc perturbations. Maxwell claimed, and 
others have confirmed, that these out-of-plane perturbations are less destabilizing than in- 
plane ones and hence our analysis, while not fully general, does get to the right answer. 
Our main reason for wishing to restrict to the planar case is that we can then work in the 
complex plane and our entire analysis can be carried out purely in complex notation, which 
makes the equations and derivations more compact, more elegant and therefore, we hope, 
more transparent. 



2. Equally-Spaced, Equal-Mass Bodies in a Circular Ring About a Massive 

Body 

Consider the multibody problem consisting of one large central body, say Saturn, having 
mass M and n small bodies, such as boulders, each of mass m orbiting the large body in 
circular orbits uniformly spaced in a ring of radius r. Indices to n — 1 will be used to denote 
the ring masses and index n will be used for Saturn. Throughout the paper we assume that 
n > 2. For the case n — 1, Lagrange proved that the system is stable for all mass ratios 
m/M. 

The purpose of this section is to show that such a ring exists as a solution to Newton's 
law of gravitation. In particular, we derive the relationship between the angular velocity 
u of the ring particles and their radius r from the central mass. We assume all bodies lie 
in a plane and therefore complex-variable notation is convenient. So, with i — ^/ —1 and 
z = X + iy, we can write the equilibrium solution for j = 0, 1, . . . , n — 1, as 

and 

Zn = 0. (2) 

By symmetry (and exploiting our assumption that n >2), force is balanced on Saturn itself. 
Now consider the ring bodies. Differentiating (1), we see that 



= -Uj'^Zj. 



(3) 
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Prom Newton's law of gravity we have that 

\ ^ ^ 



-GM P \ + y Cm , ' (4) 



Equations (3) and (4) allow us to determine a;, which is our first order of business. By 
symmetry it suffices to consider j — Q. It is easy to check that 

Zk-ZQ = re^'^*e"^'=/"2i sin(7rA;/n) (5) 

and hence that 

\zk ~ Zo\ — 2r sin(7r/c/n). (6) 
Substituting (5) and (6) into (4) and equating this with (3), we see that 

GM Gm ze^^^/" 

^ 4'^^ sin^(7r/c/n) 

CM Gm^ 1 .Gm ^ cos(7rA;/n) 



El .Gm cosyKk/n) , . 

sin ^'iri' /n zLt3 Z-^ a\-r\'^ ( I n\ 



J.3 4^3 sin(7rA;/n) 4r^ sin^(7rA;/nj 

It is easy to check that the summation in the imaginary part on the right vanishes. Hence, 

2 GM Gm^ 

where 

1 1 

^- = 4S^sin(7rfc/n)- 

With this choice of a;, the trajectories given by (1) and (2) satisfy Newton's law of gravitation. 



3. First-Order Stability 

In order to carry out a stability analysis, we need to counter-rotate the system so that 
all bodies remain at rest. We then perturb the system slightly and analyze the result. 

A counter-rotated system would be given by 

In such a rotating frame of reference, each body remains fixed at its initial point. It turns 
out to be better to rotate the different bodies different amounts so that every ring body is 
repositioned to lie on the a;-axis. In other words, for j = 0, . . . , n — 1, n, we define 

wj = Uj + ivj = e-^^'^*+2'^^/") Zj. (11) 
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The advantage of repositioning every ring body to the positive real axis is that perturbations 
in the real part for any ring body represent radial perturbations whereas perturbations in 
the imaginary part represent azimuthal perturbations. A simple counter-rotation does not 
provide such a clear distinction between the two types of perturbations (and the associated 
stability matrix fails to have the circulant property that is crucial to all later analysis). 

Differentiating (11) twice, we get 
Prom Newton's law of gravity, we see that 

where 



Wj — to Wj 



2iu;Wj + J2Gmkjf^, (13) 



Jm, forA; = 0, l,...,n-l, 
'^'=\m, iork = n, 

Cfcj = e'^^-^Wk - Wj (15) 

and 

9k = 2nk/n. (16) 

Let Swj{t) denote variations about the fixed point given by 

fr for, = 0,l,...,n-l, 
[ 0, for J = n. 

We compute a linear approximation to the differential equation describing the evolution of 
such a perturbation. Applying the quotient, chain, and product rules as needed, we get 

= A, - 2i.5w, - l^Gm J^^'^''^^;'! (18) 

where 
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The next step is to use (15) to re-express the ^k/s in terms of the WkS and the Wj^s and 
then to substitute in the particular solution given by (17). Consider the case where j < n. 
In this case we have 

r(e*^'=-j — 1), for /c < n, 
— r, for A; = n 



and therefore 



2rsin(|6'fe_j-|/2), ior k < n, 
r, for k — n. 



Substituting these into (18) and simplifying, we get 

CM, , o ,« . X GM 
^(e-^^5«;„ + 3e^^^5^„) + ^ 



4. Choice of Coordinate System 

Without loss of generality, we can choose our coordinate system so that the center of 
mass remains fixed at the origin. Having done that, the perturbations 6wn and 6tVn can be 
computed explicitly in terms of the other perturbations. Indeed, conservation of momentum 
implies that 

rn^Szk + MSzn = 0. 

kj^n 

Hence, 

From the definition (11) of the WkS in terms of the ZkS, it then follows that 

^-'''^^n--^Y.^''''-^^Wk. 

ky^n 

Making this substitution for e~^^^Swn and an analogous substitution for e^^^Swn in (19), we 
see that 

5wj = u^Swj - 2iuj5wj + y^T^ {e'^''-'5wk + Se"'^''-^ 5wk) + ^^{Swj + 3(5%) 

k^n 

Gm 1 e'-^^-i^Wk — Swj — ?>e^^*"'{e~''^*'-^5wk — Swj) , , 



-n3(|^,_,|/2) 
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5. Circulant Matrix 



Switching to matrix notation, let Wj denote a shorthand for the column vector [ Wj Wj ] . 
In this notation, we see that (20) together with its conjugates can be written as 



d 
di 



SWo 




5Wi 
















_ SWn-l _ 





D 

Nn-l 



D 



Nn-l 
Nn-2 



Ni N2 ■■■ D 
where D,Q,, and the A^^'s are 2x2 complex matrices given by 





5Wo 























D 



1 1 
1 1 



+ 



Gm 
2^ 



1 - /„ + Jn/2 3 - 3 J„/2 

3 - 3 Jn/2 1-In + Jn/2 



(21) 



Gm 
2^ 



e'^^ (1 - Jfc,„/2) 3e-'^'= + 3 Jfc,„/2 
3e^^'=+3Jfc,„/2 e-'''^{l-Jk,n/2) 



and where 



and 



2iu! 



-1 
1 



n— 1 



k=l 



Jk,n 



1 



4sin {7r\k\/n) 
1 

4sin^ {n\k\/ri) 



^ ("-l)/2 ^ ^ 

Y,Ik,n ~ E ^ ~ — nlog(n/2) 

fe=i 



(22) 



n— 1 oc ^ 



n 



k=l 



27r3 ^ 27r3 

k=l 



C(3) = 0.01938 



(23) 
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Here, the symbol is used to indicate asymptotic agreement. That is, a„ ~ 6„ means that 
O'n/bn ^ 1 as n — > oo and ({3) denotes the value of the Riemann zeta function at 3. This 
constant is known as Apery's constant (see, e.g., Arfken (1985)). 

Finally, note that in deriving (21) from (20) we have made repeated use of the following 
identity 



^sin^ 1^,1/2 



4J„, - 8/„ 



Let A denote the matrix in (21). We need to find the eigenvalues of A and derive 
necessary and sufficient conditions under which none of them have a positive real part. At 
this point we could resort to numerical computation to bracket a threshold for stability by 
doing a binary search to find the largest value oi m/M for which none of the eigenvalues 
have positive real part. We did such a search for some values of n. The results are shown in 
Table 1. 

The eigenvalues are complex numbers for which there are nontrivial solutions to 

/ 



(24) 



D 

Nn-l 



D 



Ni N2 



D 





SWo 




5Wo 




5Wi 








SWo 


= A 


SWn-1 

6W0 




5Wi 




5Wi 




_ SWn-l _ 







The first group of equations (above the line) can be used to ehminate the "derivative" 
variables from the second set. That is. 



5Wo 




5Wo 








= A 


. SWn-1 . 







and therefore 



D Ni ■ 






SWo 








SWo 




SWo 


Nn-l D ■ 






SWi 


+A 


n 




SWi 


= \' 


SWi 


Ni N2 ■ 


• D 












. SWn-1 . 




. SWn-1 _ 



(25) 
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The matrix on the left-hand side is called a block circulant matrix. Much is known about 
such matrices. In particular, it is easy to find the eigenvectors of such matrices. For general 
properties of block circulant matrices, see Tee (2005). 

Let p denote an n-th root of unity (i.e., p — e^'^'-'/" for some j = 0, 1, . . . , n — 1) and let 
^ be an arbitrary complex 2-vector. We look for solutions of the form 



6Wo 






5Wi 




Pi 









Substituting such a guess into (25), we see that each of the n rows reduce to one and the 
same thing 

{D + pN, + ... + p""^iV„_i) ^ + m = A^e. 
There are nontrivial solutions to this 2x2 system if and only if 

det(D + pA^i + . . . + p"~^Nn-i + = 0. 

For each root of unity, p, there are four values of A that solve this equation (counting 
multipli cites) . That makes a total of in eigenvalues and therefore provides all eigenvalues 
for the full system (24). 



6. Explicit Expression for Yl!k=iP 

In order to compute the eigenvalues, it is essential that we compute Yl^^P^^k as 
explicitly as possible. To this end, we note the following reduction and new definition: 



n— 1 ^ n— 1 



k=l 



1 Y^_cos(j6'fc) 



4ttsin^(^fc/2) 
Jj,n- (26) 



Similarly, 



n-l 



J2p'^'''Jk,n^Jj+l,n (27) 



k=l 



and 



Also we compute 



n-l 



n-1 



^ ] p e ^ Jk,n — Jj—\,n- 



k=l 



n-l 



n-l 



{j+i)efc 



fe=i 



ifc=i 



ifc=i 



n — 1 for j = n — 1 
— 1 otherwise 



and 



n-l 



k=l 



n — 1 for j = 1 
— 1 otherwise . 



Substituting the definition of A^fc into Yl^=iP''^k and making use of (26)-(30), we get 



n-l 



fe=l 



Gm 
2^ 



-1 + nSj=n-i - lJj+i,n -3 + 3nSj=i + \ jj^n 
—3 + 2>n5j=n-i + \ Jj,n ~1 + 'ndj=i — | Jj_i^„ 



where (5^=^ denotes the Kronecker delta (i.e., one when j — k and zero otherwise). 

7. Solving det {D + Y2Z\p^Nk + Afi - = 0. 
Assembling the results from the previous sections, we see that 

n-l 

D + Y P^^k + AO - 



k=l 



+ - - 2iu}X - A^ 



3, ,2 3^2 
2^ - 2°^j 



3, ,2 3^2 



|a;2+iQ,2_^_/32^2iu;A-A2 



Cm 
2^ 



3n5j=n^i n5j=i 



where and /3'^ are shorthands for the expressions 



a. 



Gm 
2^ 



and 
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and, as a reminder, 7„ and J„ are defined by (22) and (23), respectively, whereas Jj,„ is 
defined by (26). 

It turns out in our subsequent analysis that the root of unity given by ^' = n/2 is the 
most critical one for stability, at least for n>7. For n — 2, ... ,6 the instability stems from 
the eigenvectors associated with j — 1 and j — n — 1. We will analyze the key cases. But 

first, wc note that the critical j = n/2 case corresponds to perturbations in which every 
other body is perturbed in the opposite direction. And, more importantly, it doesn't matter 
what the direction of the perturbation is. That is, if body is advanced azimuthally, then all 
of the even-numbered bodies are advanced azimuthally and all of the odd-numbered bodies 
are retarded by the same amount. Similarly, if body is pushed outward radially, then all 
of the even-numbered bodies are also pushed outward whereas the odd-numbered bodies are 
pull inward. Azimuthal and radial perturbations contribute equally to instability. 



7.1. The Case Where n Is Arbitreiry And j Is Neither 1 Nor n — 1. 

Assuming that j is neither 1 nor n — 1, we see that 

det (^D + i^P^'^fc + Xn- A^/j = Qcu^ + -(3^- 2iujX - A^^ 



3 o 1 
2' 



Expanding out the products on the right-hand side in (33), we get that 



det Ld + J^p'^^k + Afi - A^J = A^ + AjX^ + iBjX + Cj = 0, (34) 



k=l 



where 



Aj = u;'-^{a'j_, + a]^,)+2p' (35) 

Bj = -a; - a|+i) (36) 

- 3c.^Q(a,^ + a|^J + ^a5-/?^) + Q(a|_, + a,V)-/?^y 

-l^«i-"l.i)^-i4 (37) 



- 12 - 



7.1.1. The Subcase Where n Is Even And j — n/2. 

In the subcase where n is even and j = n/2, it is easy to see by symmetry that an/2+1 = 
<^n/2-i- To emphasize the equahty, we will denote this common value by an/2±i- Equations 
(35) to (37) simplify significantly. The result is 



n-l 



det LD + J^P'^^k + Xn-X^l\ = + {u^ - a^/2±i + 2p^)X^ 



k=l 



1 3 
+3lj^ ( -an/2±i + -an/2 - /^^ 



1 Y 9 

2"n/2±l-/^M -4"n/2- (38) 

For a moment, let us write this biquadratic polynomial (in A) in a simple generic form and 
equate it to zero 

A^ + AA^ + C = 0. 
The quadratic formula then tells us that 

2 ■ 

To get the eigenvalues, we need to take square roots one more time. The only way for 
the resulting eigenvalues not to have positive real part is for A^ to be real and nonpositive. 
Necessary and sufficient conditions for this are that 

(39) 
(40) 
(41) 

It turns out that the third condition implies the first two (we leave verification of this fact 
to the reader). In terms of computable quantities, this third condition can be written, after 
simplification, as 

+ i-8al/,^, - 18al/, + 16/3')uj' + 9a^/2 > 0. 
Again we use the quadratic formula to find that 



A 


> 





C 


> 





AC 


> 


0. 



or 



< 4«n/2±l + 9«n/2 ' ^P' ' ^(4«n/2±l + 9<2 ' ' 



4 

n/2- 



4 

n/2- 
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It is the greater-than constraint that is relevant and so we take the positive root. Finally, 
we recall that 



2 

= — ^ H 5-4 

2 _ Gm ~ 

— l2^\Jn- Jn/2,n) 

and so the inequality on o;^ reduces to 

M . 9 _ 

> 2( J„ — J„/2±i n) + 7:(<4 ~ ■4/2,n) ~ 5/„ 

+ ^2(J„ - Jn/2±l,n) + 2('^" ~ 4/2,71) " - - - Jn/2,n^ ■ (42) 

The second column in Table 1 shows thresholds computed using this inequality It is 
clear that for even values of n greater than 7, this threshold matches the numerically derived 
threshold shown in the first column in the table. This suggests that inequalities analogous 
to (42) derived for j ^ n/2 are less restrictive than (42). The proof of this statement is 
obviously more complicated than the j = n/2 case because the general case includes a linear 
term {Bj ^ 0) which vanishes in the j = n/2 case. The linear term makes it impossible simply 
to use the quadratic formula and therefore any analysis involves a more general analysis of 
a quartic equation. Scheeres and Vinh (1991) analyzed the general case. Although their 
notations are different, the fundamental quantities are the same and so their analysis is valid 
here as well. Rather than repeating their complete analysis, we simply outline the basic 
steps in the next subsection. 



GM Gm 



7.1.2. The Subcase Where j^n/2. 

Let aj±i denote the average of aj+i and aj^i. 

otj±i = (ftj+i + Q;j-i)/2. 

If we were to assume, incorrectly, for the moment that terms involving the difference ctj+i — 
otj-i were not present in (35)-(37), then an analysis analogous to that given in the previous 
subsection would give us the following inequality: 
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Next, one uses the fact that o:^±i is unimodal as a function of j taking its maximum value at 
j = n/2. Hence, the inequahty associated with j = n/2 is the strictest of these inequahties. 
Finally, the difference terms are treated as small perturbations to this simple case and a 
homotopy analysis shows that the j — n/2 case remains the strictest case even as the 
difference terms are fed in. 



7.1.3. Large n 

When n is large, J„/2±i,n ~ Jn/2,n and In- Furthermore, 



J, 



n/2,n 



2f^sm^(kTT/n) ~ 2%^ f-^ 

k=l ^ ' ' fc=l 



3 ^3 - 1 3 1 1 3 



1 _ 3 1 ^ 



Hence, (42) reduces to 
or, equivalently. 



4 27r3 A;3 4 2 sin3(A;7r/n) 4 " 
fc=i fc=i ^ ' ' 



— > - 13 + 2^ J„, 
m 8 

m < ^^-^ 2.299M/n^ (43) 

- 1(13 + 2^) J„ ^ ' ^ ^ 

which is precisely the answer Maxwell obtained 150 years ago. Of course, we have assumed 
here that n is even. For the odd case, as n — > oo, — CKj+i] — > so that the odd quartic 
equation for j — {n — l)/2 reduces to the the even equation for j — n/2 giving the same 
stability criteria as the even particle case. This can be seen in our simultions as well. The 
case n = 100 and n = 101 give the same threshold to several significant figures. 



7.2. The Case Where j = 1 Or j ^n-1. 

By symmetry, these two cases are the same. Hence, we consider only j — 1. Again after 
some manipulation in which we exploit the fact that ccq = 0, we arrive at 

det (d + J^p'^^k + An - A^/ ) = A^ + AjX'^ + iBjX + Cj = 0, (44) 



k=l 



where 



Ai = cu^ -^{n-f + al) +2p^ (45) 
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Bi = -uj{n-f-al) (46) 



and 



-^{wy-alf -^al{al-n-f), (47) 



(48) 



Note that the coefficient Bi is imaginary whereas the other three coefficients are real. 
This suggests making the substitution jj, = iX. In terms of n, (44) becomes a quartic equation 
with all real coefficients: 

/x^ - Aiii^ + Bi^ + Ci = 0. (49) 

This equation either has four real roots or not. If it does, then the corresponding values 
for A are purely imaginary and the system could be stable. If, on the other hand, there are 
two or fewer real roots, then at least one pair of roots to (49) form a conjugate pair and 
therefore the corresponding pair of values for A will be such that one has positive real part 
and the other negative real part. Hence, in that case, the system is demonstrably unstable. 
Simple numerical investigation reveals that this is precisely what happens when 2 < n < 6 
regardless of the mass ratio M/m. 

To see why, let us consider just the case when M/m is very large and hence the ratio 

r = m/M 

is very close to zero. In this asymptotic regime, 

Ai = aM 
Bi = bVMm 
Ci = cMm, 

where a > and the sign of c is the same as the sign of |q;| + |Q;f — — ^n'y. Substituting 
rM for m and making the change of variables defined hy u — /i/y/M, we get 

i/^ — au^ + hrv + cr = 0. 

For r = 0, this equation reduces to i/^ — af'^ — which has three real roots, a positive 
one, a negative one, and a root of multiplicity two at z/ = 0. By continuity, for r small but 
nonzero, the quartic still has a positive root and a negative root but the double root at = 
can either disappear or split into a pair of real roots. Since this bifurcation takes place in 
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a neighborhood of the origin, the quartic term can be ignored and the equation in a tiny 
neighborhood of zero reduces to a quadratic equation: 

—av^ + hrv + cr = 0. 

This equation has two real roots if and only if its discriminant is nonnegative: 

r(r6^ + 4ac) > 0. 

Hence, if c is negative there will not be a full set of real roots for r very small and hence the 
ring system will be unstable in that case. In other words, the system will be unstable if 

+ \(A - - \ni < 0. 

This equation reduces to 

n-l ^ ^ 

V ^-rrr - n - - cot (^) < 0. (50) 

f^sin(^) 2 \2nJ ^ ^ 

k=l \ n / 

It is easy to check that the expression is negative for n = 2, . . . , 6 and positive for n > 7. 
Therefore, we have proved that ring systems are unstable for n = 2, . . . , 6 at least when m is 
very small relative to M. We have not proved the result for larger values of m but it seems 
that such a case should be even more unstable, which is certainly verified by our simulator. 



8. Ring Density 

Suppose that the linear density of the boulders is A. That is, A is the ratio of the 
diameter of one boulder to the separation between the centers of two adjacent boulders. 
Then the diameter of a single boulder is X{2nr/n). Hence, the volume of a single boulder is 
(47r/3)(A7rr/n)^. Let 5 denote the density of a boulder. Then the mass of a single boulder is 
(47r/3)(A7rr/n)^(5. If we assume that the density of Earth is about 8 times that of a boulder 
(Earth's density is 5.5 and Saturn's moons have a density of about 0.7 being composed of 
porous water-ice), then we have 

"8(47r/3)r|' 

where niE denotes the mass of Earth and denotes its radius. Combining all of these 
factors and assuming the central mass is equal to Saturn's mass and the ring's radius is 
about the radius of the Cassini division (120, 000km), we see that the upper bound on the 
linear density of boulders is 

^-l,^25.65)(0.01938)j V ' ^'^^^^ 
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In other words, the hnear density cannot exceed 22% otherwise the ring will be unstable. Of 
course, this is for a one dimensional circular ring of ice boulders. Analysis of a two dimen- 
sional annulus or the full three dimensional case is naturally more complicated. Nonetheless 
the 22% linear density figure matches surprisingly well with the measured optical density 
which hovers around 0.05 to 2.5. 



9. Numerical Results 

We have computed stability thresholds three different ways for various finite n. 

First, we numerically solved for all eigenvalues of the Ax x An matrix in (24) and did 
a binary search to locate the smallest mass ratio M/ m for which no eigenvalue has positive 
real part. We then tranlated this threshold into a value of 7 for the threshold expressed as 

m < -iM/n^ 

and tabulated those results in the column labeled Numerical in Table 1. 

Secondly, for even values of n we used equation (42) to derive 7 threshold values. These 
values are reported in the second column of thresholds in Table 1. Note that for even values 
of n larger than 7, these results agree with those obtained numerically. 

Lastly, and perhaps most interestingly, the third column of results in the table are 
stability thresholds that were estimated using a simulator (Vandcrbci (2005)) based on a 
leap-frog integrator (Saha and Tremaine (1994); Hut et al. (1995)). In this column, two 
values are given. For the larger value instability has been decisively observed. However, 
verifying stability is more challenging since one should in principle run the simulator forever. 
Rather than waiting that long, we use the rule of thumb that if the system appears intact for 
a period of time ten times greater than the time it took to demonstrate instability, then we 
deem the system stable at that mass. This is how the lower bounds in the table were obtained. 
It is our belief that a good (simplectic) simulator provides the most convincing method to 
discriminate between stable and unstable orbits at least when the number of bodies remains 
relatively small, say up to a few hundred. Unstable orbits reveal themselves quickly as the 
initial inaccuracy of double precision arithmetic quickly cascades into something dramatic if 
the system is unstable. If, on the other hand, the system is stable then the initial imprecisions 
simply result in an orbit that it close to but not identical to the intended orbit. The situation 
does not decay. Any reader who has never experimented with a good simplectic integrator 
is strongly encouraged to experiment with the Java applet posted at 

http: / /www.princeton.edu/ rvdb/JAVA/astro/galaxy/StableRings.html 
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as hands on experience can be very convincing. 

Of course, the amazing thing about the simulator results is that they match the numer- 
ical results in the first column. The thresholds determined by linear stability analysis only 
tell us definitively that for m larger than the threshold, the system is necessarily unstable. 
But, for m smaller than the threshold, the mathematical/numerical analysis says nothing 
since in those cases, the eigenvalues are all purely imaginary. Yet, simulation confirms that 
the thresholds we have derived are truely necessary and sufficient conditions for stability. 

As shown in Section 7.2 for 2 < n < 6, the system is unstable. The simulator verifies 
this. In these cases, there is lots of room to roam before one body catches up to another. 
Even the tiniest masses are unstable. The case n = 2 is especially interesting. For this case, 
the two small bodies are at opposite sides of a common orbit. Essentially they are in L3 
position with respect to each other. For the restricted 3-body problem (where one body has 
mass zero) it is well-known that L3 is unstable no matter what the mass ratio. Our results 
bear this out. 
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n 


Threshold Value for 7 

Numerical Eq. (42) Simulator 
(even n only) 


2 


0.000 


4 


[0.0, 0.007] 


6 


0.000 


2.487 


[0.0, 0.025] 


7 


2.452 




[2.45, 2.46] 


8 


2.412 


2.4121 


[2.41, 2.42] 


10 


2.375 


2.3753 


[2.37, 2.38] 


36 


2.306 


2.3066 


[2.30, 2.35] 


100 


2.300 


2.2999 


[2.30, 2.31] 


101 


2.300 




[2.30, 2.31] 



Table 1: Estimates of the stability threshold (i.e., 7 in an inequality of the type m < 'jM/n^). 
The first column contains numerically derived obtained by a brute-force computation of 
the eigenvalues together with a simple binary search to find the first point is at which an 
eigenvalue takes on a positive real part. The second column gives thresholds computed 
using (42). The column of simulator values corresponds to results from running a leap-frog 
integrator and noting the smallest value of 7 for which instability is clearly demonstrated. 
This is the larger of the pair of values shown. The smaller value is a nearby value for which 
the simulator was run ten times longer with no overt indication of instability. 
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